##### ---- Figure 1c and 1d Replication ---- #####

rm(list = ls())
# Load required libraries
library(tidyverse)
library(cowplot)
library(robustbase)
library(matrixStats)


# Set working directory
setwd("")

# Set seed
set.seed(1381) # Same as Wilkins

# Parameters for the DGP
phi2 <- seq(0, 0.5, 0.1)
beta <- 0.5
alpha <- 0.75
rho <- 0.95
n <- 100
sims <- 1000 # No. of MC sims

# List to store results
mat_coef <- list()
mat_pvals <- list()
mat_lb <- list()
mat_ub <- list()
mat_se <- list()

# Monte carlo analysis, figure 1c and 1d
for(j in 1:length(phi2)){
  phi <- phi2[j]
  
  # Matrices to store results
  mat_coef[[j]] <- matrix(data = NA, nrow = sims, ncol = 4)
  colnames(mat_coef[[j]]) <- c("eq4", "lgdv", "lgdv2", "reg")
  mat_se[[j]] <- matrix(data = NA, nrow = sims, ncol = 4)
  colnames(mat_se[[j]]) <- c("eq4", "lgdv", "lgdv2", "reg")
  mat_pvals[[j]] <- matrix(data = NA, nrow = sims, ncol = 4)
  colnames(mat_pvals[[j]]) <- c("eq4", "lgdv", "lgdv2", "reg")
  mat_lb[[j]] <- matrix(data = NA, nrow = sims, ncol = 4)
  colnames(mat_lb[[j]]) <- c("eq4", "lgdv", "lgdv2", "reg")
  mat_ub[[j]] <- matrix(data = NA, nrow = sims, ncol = 4)
  colnames(mat_ub[[j]]) <- c("eq4", "lgdv", "lgdv2", "reg")
  
  for(i in 1:sims) {
    # Generate autoregressive x
    x <- numeric(n+2)
    x[1] <- rnorm(1, 0, 1)
    for(t in 2:length(x)){
      x[t] <- rho*x[t-1] + rnorm(1, 0, 1)
    }
    
    # Generate u (serially correlated errors)
    u <- numeric(n+2)
    u[1] <- rnorm(1, 0, 1)
    for(t in 2:length(u)){
      u[t] <- phi*u[t-1] + rnorm(1, 0, 1)
    }
    
    # Generate y
    y <- numeric(n+2)
    yinit <- rnorm(1, 0, 1)
    
    for(t in 1:length(x)){
      if(t==1) y[t] <- alpha*yinit + beta*x[t] + u[t]
      if(t > 1) y[t] <- alpha*y[t-1] + beta*x[t] + u[t]
    }
    
    ###Prepare time series and lagged variable data
    ynl <- y[-(1:2)]
    ylag1 <- y[-c(1, length(y))]
    ylag2 <- y[-(length(y):(length(y)-1)) ]
    
    xn <- x[-1]
    xlag <- x[-length(x)]
    
    xn <- xn[-1]
    xlag <- xlag[-1]
    
    eq4 <- lm(ynl ~ 0 + ylag1 + ylag2 + xn + xlag)
    lgdv <- lm(ynl ~ 0 + ylag1 + xn)
    lgdv2 <- lm(ynl ~ 0 + ylag1 + ylag2 + xn)
    reg <- lm(ynl ~ 0 + xn)
    
    # EQ4
    mat_coef[[j]][i, 1] <- summary(eq4)$coefficients[3, 1]
    mat_se[[j]][i, 1] <- summary(eq4)$coefficients[3, 2]
    mat_pvals[[j]][i, 1] <- summary(eq4)$coefficients[3, 4]
    mat_lb[[j]][i, 1] <- confint(eq4)[3, 1]
    mat_ub[[j]][i, 1] <- confint(eq4)[3, 2]
    
    # LGDV
    mat_coef[[j]][i, 2] <- summary(lgdv)$coefficients[2, 1]
    mat_se[[j]][i, 2] <- summary(lgdv)$coefficients[2, 2]
    mat_pvals[[j]][i, 2] <- summary(lgdv)$coefficients[2, 4]
    mat_lb[[j]][i, 2] <- confint(lgdv)[2, 1]
    mat_ub[[j]][i, 2] <- confint(lgdv)[2, 2]
    
    # LGDV2
    mat_coef[[j]][i, 3] <- summary(lgdv2)$coefficients[3, 1]
    mat_se[[j]][i, 3] <- summary(lgdv2)$coefficients[3, 2]
    mat_pvals[[j]][i, 3] <- summary(lgdv2)$coefficients[3, 4]
    mat_lb[[j]][i, 3] <- confint(lgdv2)[3, 1]
    mat_ub[[j]][i, 3] <- confint(lgdv2)[3, 2]
    
    # REG
    mat_coef[[j]][i, 4] <- summary(reg)$coefficients[1, 1]
    mat_se[[j]][i, 4] <- summary(reg)$coefficients[1, 2]
    mat_pvals[[j]][i, 4] <- summary(reg)$coefficients[1, 4]
    mat_lb[[j]][i, 4] <- confint(reg)[1, 1]
    mat_ub[[j]][i, 4] <- confint(reg)[1, 2]
  }
}


# Perecent Bias
bias <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
  bias[i, ] <- (colMeans(mat_coef[[i]] - beta))/beta * 100
}
bias <- cbind(phi2, bias)
colnames(bias)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
bias <- as.tibble(bias) 
bias <- bias %>%
  select(-REG) %>%
  gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "bias")

percentbias <- ggplot(bias, aes(x = phi, y = bias, color = Model)) +
  geom_line(aes(linetype = Model), size = 1) + 
  ylim(-50, 50) +
  ylab("") + 
  xlab(expression(phi)) +
  ggtitle("Percent Bias") + 
  theme_minimal() +
  theme(legend.position = "none") +
  scale_linetype_manual(values=c("solid", "longdash", "dotdash"))

# RMSE
rmse <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
  rmse[i, ] <- sqrt(colMeans((mat_coef[[i]] - beta)^2))
}
rmse <- cbind(phi2, rmse)
colnames(rmse)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
rmse <- as.tibble(rmse)
rmse <- rmse %>%
  select(-REG) %>%
  gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "rmse")

rootmean <- ggplot(rmse, aes(x = phi, y = rmse, color = Model)) +
  geom_line(aes(linetype = Model), size = 1) + 
  ylab("") + 
  ylim(0, 0.25) +
  xlab(expression(phi)) +
  ggtitle("RMSE") + 
  theme_minimal() +
  theme(legend.position = "none") +
  scale_linetype_manual(values=c("solid", "longdash", "dotdash"))

# Power 
power <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
  power[i, ] <- colMeans(ifelse(mat_pvals[[i]] <= 0.05, 1, 0))
}
power <- cbind(phi2, power)
colnames(power)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
power <- as.tibble(power)
power <- power %>%
  select(-REG) %>%
  gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "power")

pow <- ggplot(power, aes(x = phi, y = power, color = Model)) +
  ylim(0.95, 1) +
  geom_line(aes(linetype = Model), size = 1) + 
  ylab("") + 
  xlab(expression(phi)) +
  ggtitle("Power") + 
  theme_minimal() +
  theme(legend.position = "none") +
  scale_linetype_manual(values=c("solid", "longdash", "dotdash"))

# Coverage 
cov <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
  cov[i, ] <- colMeans(ifelse(mat_lb[[i]] <= beta & mat_ub[[i]] >= beta, 1, 0))
}
cov <- cbind(phi2, cov)
colnames(cov)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
cov <- as.tibble(cov)
cov <- cov %>%
  select(-REG) %>%
  gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "coverage")

cover <- ggplot(cov, aes(x = phi, y = coverage, color = Model)) +
  geom_line(aes(linetype = Model), size = 1) + 
  ylim(0, 1) +
  xlab(expression(phi)) +
  ylab("") + 
  ggtitle("Coverage") + 
  theme_minimal() +
  theme(legend.position = "none") +
  scale_linetype_manual(values=c("solid", "longdash", "dotdash"))

# Standard Deviation
sd <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
  sd[i, ] <- colSds(mat_coef[[i]])
}
sd <- cbind(phi2, sd)
colnames(sd)<- c("phi","EQ4", "LGDV", "LGDV2", "REG")
sd <- as.tibble(sd)
sd_wide <- sd
sd <- sd %>%
  select(-REG) %>%
  gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "sd")

stddev <- ggplot(sd, aes(x = phi, y = sd, color = Model)) +
  geom_line(aes(linetype = Model), size = 1) + 
  ylim(0, 0.2) +
  xlab(expression(phi)) +
  ylab("") + 
  ggtitle("Standard Deviation") + 
  theme_minimal() +
  theme(legend.position = "none") +
  scale_linetype_manual(values=c("solid", "longdash", "dotdash"))

# Optimism
mean_se <- matrix(NA, 6, 4)
for(i in 1:length(phi2)){
  mean_se[i, ] <- colMeans(mat_se[[i]])
}
mean_se <- cbind(phi2, mean_se)
colnames(mean_se)<- c("phi","se_EQ4", "se_LGDV", "se_LGDV2", "se_REG")
mean_se <- as.tibble(mean_se) 

opt <- left_join(sd_wide, mean_se, by = "phi")
opt <- opt %>% mutate(overconf_EQ4 = EQ4/se_EQ4,
                      overconf_LGDV = LGDV/se_LGDV,
                      overconf_LGDV2 = LGDV2/se_LGDV2,
                      overconf_REG = REG/se_REG)

opt <- opt %>% select(phi, overconf_EQ4, overconf_LGDV, overconf_LGDV2, overconf_REG)
colnames(opt) <- c("phi","EQ4", "LGDV", "LGDV2", "REG")
opt <- opt %>%
  select(-REG) %>%
  gather(`EQ4`, `LGDV`, `LGDV2`, key = "Model", value = "overconf")


optimism <- ggplot(opt, aes(x = phi, y = overconf, color = Model)) +
  geom_line(aes(linetype = Model), size = 1) + 
  ylim(0, 2) +
  xlab(expression(phi)) +
  ylab("") + 
  ggtitle("Overconfidence") + 
  theme_minimal() +
  theme(legend.position = "none") +
  scale_linetype_manual(values=c("solid", "longdash", "dotdash"))

# Combine plots
plot_grid(rootmean, cover, pow, percentbias, stddev,  optimism, nrow = 3)
ggsave("WilkinsRep.pdf", width = 8.5, height = 11, units = "in")


